      !<@file 
      !>@brief PDMDFDT1: Problema Direto, Método das Diferenças Finitas, Domínio do Tempo, 1: 2a ordem no tempo e 4a ordem no espaço.
      !>
      !> Baseado na implementação feita na disciplina de métodos sísmicos
      subroutine PDMDFDT1(nttotal,ntps,St,nx,nz,delta,deltat,c,rho,&
		    & isourcefixa,jsourcefixa,Nrec,ireceiv,jreceiv,&
		    & P1,P2,P3,Aaux,Caux,plotsnap,ntsnaps,PtrecR)

      implicit none


      integer :: nttotal !< número total de instantes de tempo 
      integer :: ntps !< amostragem, número de instantes de tempo por segundo
      double precision :: St(nttotal) !< fonte no domínio do tempo
      integer :: nx !< número de nós no grid de diferenças finitas na direção x
      integer :: nz !< número de nós no grid de diferenças finitas na direção z
      double precision :: delta !< Distância entre os nós da malha de diferenças finitas nas direções x e z
      double precision :: deltat !< intervalo de tempo da discretização temporal
      double precision :: rho(nx,nz) !< Densidade (massa específica)
      double precision :: c(nx,nz) !< Velocidade de propagação
      integer :: isourcefixa !< Número do nó na direção x da fonte referente a esta simulação
      integer :: jsourcefixa !< Número do nó na direção z da fonte referente a esta simulação
      integer :: Nrec !< Número de receptores referentes à fonte em simulação
      integer :: ireceiv(Nrec) !< Número do nó na direção x dos receptores.
      integer :: jreceiv(Nrec) !< Número do nó na direção z dos receptores.
      integer :: plotsnap !< plotsnap = 1 salva o campo em todo domínio em bin para futuro snapshot
      integer :: ntsnaps !< Número de instantes de tempo gravados para os snapshots.
      real :: PtrecR(nttotal,Nrec) !< campo de onda em cada instante de tempo em cada receptor



      integer :: itsnap !< Número do próximo instante de tempo que será gravado para snapshot



      integer :: i
      integer ::j




      integer :: it !< identifica o passo de tempo
      integer :: kr !< identifica o número do receptor
      integer :: ierr !< erro de alocação

      double precision :: P1(nx,nz),P2(nx,nz),P3(nx,nz)
      double precision :: Caux(nx,nz),Aaux(nx,nz)

      integer :: iplotsnap!< contador de plotagens em bin para os snapshots


      ! inicialização do campo antes do loop temporal
      do i=1,Nx
	  do j=1,Nz
	      P1(i,j) = 0.d0
	      P2(i,j) = 0.d0
	      P3(i,j) = 0.d0
	  enddo
      enddo 

      iplotsnap = 1 ! inicializando contador de plotagens em bin para snapshots
      itsnap = 1 ! próximo instante que terá o campo plotado
      



      do it=1,nttotal! Loop temporal

	  P3 = 0.0d0

	  ! fonte
! 	      P2(isourcefixa,jsourcefixa) = -rho(isourcefixa,jsourcefixa)*St(it-1)
	      P3(isourcefixa,jsourcefixa) = - (c(isourcefixa,jsourcefixa) * deltat)**2*St(it)/delta**2
! 	      P3(isourcefixa,jsourcefixa) = St(it)

	  !calculo do campo no interior do modelo
	  do j=3,nz-2
	      do i=3,nx-2
		  P3(i,j) = P3(i,j) + Caux(i,j)*(     P2(i+2,j) + P2(i-2,j) + P2(i,j+2) + P2(i,j-2) &
					  &  -16.d0* (P2(i+1,j) + P2(i-1,j) + P2(i,j+1) + P2(i,j-1))&
					  & + 60.d0* P2(i,j)  ) + 2.d0*P2(i,j) - P1(i,j)
	      enddo
	  enddo


	  ! cond de bordas
	  do i=3,nx-2

	      do j=1,2
! 		     p3(i,j)=p2(i,j)+Aaux(i,j)*(p2(i,j+1)-p2(i,j))
		     p3(i,j)=0.0d0
	      enddo

	      do j= nz-1,nz
		  p3(i,j)= p2(i,j)-Aaux(i,j)*(p2(i,j)-p2(i,j-1))
	      enddo
	  enddo

	  do j=3,nz-2
	      do i= 1,2
		  p3(i,j)=p2(i,j)+Aaux(i,j)*(p2(i+1,j)-p2(i,j))
	      enddo

	      do i= nx-1,nx
		  p3(i,j)=p2(i,j)-Aaux(i,j)*(p2(i,j)-p2(i-1,j))
	      enddo
	  enddo




	  ! salvando sismogramas
	  do kr=1,Nrec
	      PtrecR(it,kr) = P3(ireceiv(kr),jreceiv(kr))
	  enddo


	  ! plotando campo de onda para snapshots, se for o caso
	  if(plotsnap.eq.1 .and. it.eq.itsnap)then
	      write(1,rec=iplotsnap)((real(P3(i,j)),j=1,Nz),i=1,Nx)

	      iplotsnap = iplotsnap  + 1 ! número da próxima plotagem
	      itsnap = itsnap + int(nttotal/(ntsnaps-1)) ! número do instante de tempo da próxima plotagem
	      if(itsnap.gt.nttotal)itsnap=nttotal
	  endif








	  ! atualizacao do campo de onda
	  P1=P2
	  P2=P3



      enddo  ! fechamento do loop temporal


     





!       do igsr=1,gruposSR !loop nos grupos de fonte/receptor
! 
! 	
! 
! 	  ! salva a resposta em todo o dominio, para todas as fontes, para a freq kw 
! 	  !(só serve para fazer snapshot)
! 	  ideltaSnap = int(Nsource(igsr)/(Nsourcesnap-1))
! 	  if (iflagF.eq.0 .and. iflagG.eq.0 .and. iflagH.eq.0)then! soh simulacao (snapshot)
! 	      do kss=1,NsourceSnap
! 		  ks=(kss-1)*ideltaSnap + 1
! 		  if(ks.gt.Nsource(igsr))ks=Nsource(igsr)
! 		  do i=1,Nx
! 		      do j=1,Nz
! 			  linha = (j-1)*Nx + i ! numero do noh i,j
! 			  P(i,j,kss,igsr) =  f(linha,ks)
! 		      enddo
! 		  enddo
! 	      enddo
! 	  endif
! 
! 
! 	  ! salva a resposta nos receptores, para todas as fontes, para a freq kw 
! 	  !(serve para fazer a inversao e p gerar sismograma)
! 	  do ks=1,Nsource(igsr)
! 	      do kr=1,Nrec(igsr)
! 		  linha = (jReceiv(kr,igsr)-1)*Nx + iReceiv(kr,igsr)  ! numero do noh do receptor i
! 		  Preceiv(kr,ks,kw,igsr) =  f(linha,ks)
! 	      enddo	
! 	  enddo
! 
! 
! 	
! 
!       enddo  !loop nos grupos de fonte/receptor



      return

      end subroutine PDMDFDT1